High-harmonic generation from arbitrarily oriented diatomic molecules 
including nuclear motion and field-free alignment 
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We present a theoretical model of high-harmonic generation from diatomic molecules. The theory 
includes effects of alignment as well as nuclear motion and is used to predict results for N2, O2, 
H2 and D2. The results show that the alignment dependence of high-harmonics is governed by the 
symmetry of the highest occupied molecular orbital and that the inclusion of the nuclear motion in 
the theoretical description generally reduces the intensity of the harmonic radiation. We compare 
our model with experimental results on N2 and O2, and obtain very good agreement. 
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I. INTRODUCTION 

High-harmonic generation (HHG) is a process 
in strong-field physics that occurs when the non- 
perturbative driving of an electron by an external field 
makes it leave its parent atom or molecule, propagate in 
the field and finally when the electron is steered back to 
its origin, recombine under the emission of a high-energy 
photon The process is coherent, so the result of the 
interaction is the production of high-frequency coherent 
radiation. The emitted signal can be used to create co- 
herent UV or XUV pulses, and this radiation may then, 
e.g., in fundamental science be utilized to probe ultrafast 
attosecond electron dynamics (see Refs. 0, 13 for recent 
reviews) . 

As compared to atoms, diatomic molecules possess ad- 
ditional degrees of freedom, namely the nuclear separa- 
tion and the orientation with respect to the laser field. 
Studies of HHG from the H J molecular ion [3- 0. El along 
with the molecules N 2 and O2 0, show that the ex- 
tra parameters of the molecules give rise to additional 
effects on the harmonic spectrum. There are two rea- 
sons why this complex behavior for molecules is of inter- 
est: (i) the extra parameters available for manipulating 
HHG from the molecules may be useful tools for con- 
trolling the intensity of the harmonics, (ii) the fact that 
the harmonic spectrum contains features reminiscent of 
the molecular structure means that the emitted signals 
provide a way to study the nuclear motion and attosec- 
ond electron dynamics of the source of the high-harmonic 
radiation HEl IHII3- 

The complete theoretical modeling of HHG is compu- 
tationally very challenging. In ab initio approaches the 
time-dependent Schrodinger equation has to be propa- 
gated in a box much larger than is the case for strong- 
field ionization studies since in the former case one has to 
make sure that the electronic wave packet does not en- 
gage with an absorbing boundary. Additional complica- 
tions come from the subsequent propagation of Maxwell's 
equations in the generating medium, and in any case 
going beyond an effective one-electron model, or other 
reduced dimensionality models, seems impossible in the 
foreseeable future. These complications underline the im- 



portance of developing accurate approximate theories. In 
particular for molecules, such models, which should be 
relatively inexpensive to evaluate, are of interest since 
investigations aiming, e.g., at an increased degree of con- 
trol of the harmonic production, require a set of eval- 
uations for different bond-lengths and orientations with 
respect to the field. In this paper we present one such 
theory. We provide a model to calculate the angular de- 
pendence of the harmonic signal from a single diatomic 
molecule of arbitrary orientation with respect to the laser 
field. In particular we include the effect of nuclear mo- 
tion (vibration), and we test the accuracy of the model 
by comparison with experiments. High-harmonic gen- 
eration from diatomic molecules has already been con- 
sidered in an extension of the Lewenstein model [l3l | in 
Ref. 0. In that work depletion from the initial state 
was considered but the nuclear motion was ignored. Re- 
cently, the strong-field approximation (SFA) for HHG in 
diatomic molecules subject to ultrashort laser pulses was 
presented, again ignoring nuclear motion, but combined 
with a careful sad dle-p oint analysis and a discussion of 
the gauge problem |14j . Here, we are concerned with the 
response to a monochromatic field, and in parallel with 
the conventions in the theory of strong-field ionization, 
we will refer to our theory as the length gauge version of 
the molecular SFA for HHG. 



The present paper is organized as follows. In Sec.lTTlwc 
present the theory for HHG from diatomic molecules in 
monochromatic linearly polarized laser fields. Section ITTT1 
contains some calculational details. In Sec. IIVI we show 
the HHG yields for N2 and O2 as a function of the orien- 
tation of the molecules with respect to the laser polariza- 
tion, and compare with experimental measurements. By 
discussing HHG signals from the isotopes H2 and D2 , we 
also see how a difference in nuclear mass affects the har- 
monic spectrum. The paper ends with a summary and a 
conclusion in Sec. Appendix lAl contains a discussion 
of the momentum space wave function of a molecule in 
the single-active-electron (SAE) approximation. 
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II. THEORY 

In our treatment of the HHG from molecules, we ex- 
tend the atomic model of Ref . , and include notions 
of nuclear motion. We consider the case, where the 
laser pulse contains several cycles such that a Floquet 
approach is suitable. As we shall see below when we 
compare with experimental results, even for probe pulses 
of 50 fs 01 an d 35 fs 0] duration this turns out to be 
an accurate approximation. Hence, using the Coulomb 
gauge and the dipole approximation, a laser field with fre- 
quency lu and period T = 2ir/u> is described by the vector 
potential A(t) — Aocos(ojt). The electric field, F(t), is 
then -dA(t)/dt, i.e., F(t) = F sin(wt) with F = wA 
[atomic units, e = % = m e = a = 1, are used throughout 
unless indicated otherwise]. 

The harmonic generation is determined by the dipole 
moment 



d(t) = (y(t)\D^\*(t)), 
where the dipole operator is defined by 

n n 



(i) 



(2) 



with e the polarization vector of the field, and where 
^(t)) is the solution to the n-electron Schrodinger equa- 
tion for the diatomic molecule 
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Hn-V, 



{n) 



*(r 1> r 2> ... > r n ,R,*)=0, (3) 



with Yi the coordinate of electron i, and R the relative 
position of the two nuclei, Hq the field-free part of the 

in) 

Hamiltonian of the molecule, and V„ the interaction of 



the molecule with the laser field, 
gauge description 



In the present length 



(n) 



F(i) 



E 1 

i=l 



(4) 



From Eq. Q it is clear, that in order to obtain the 
dipole moment, we need a solution of Eq. @. The state 
IvE'(i)) is determined by using the description of the laser- 
molecule system given in Ref. [Tt| . Here nuclear mo- 
tion is modelled within the Born-Oppenheimer (BO) ap- 
proximation, and, assuming negligible electron-electron 
correlation, the electrons are put into orbitals obtained 
from Hartree-Fock (HF) calculations. Consequently, the 
interaction of the molecule can be treated in the SAE 
approximation, where the active electron is the electron 
in the highest occupied molecular orbital (HOMO). We 
write ^(t)) as a superposition of the field- free state of 
the molecule and final states, where the active electron is 
excited to the continuum. The expansion coefficients are 
obtained under the assumption that there is only limited 
ionization, which again requires that the laser intensity 
is below saturation (the regime of relevance in the HHG 
process). Inserting the expanded state into Eq. Q, we 
arrive at the following approximation to the dipole mo- 
ment 



doit) 



dt'Y, J dq<*o(t)|^ n) |*, 1 , q (*)><W* , )|vi n) (f) 1*0(0) 



(5) 



r 



Here 



det [0i(ri;i? o )02(r2;-Ro) 



(r„;-R )] 



x Xu (R) e^cp{-iE t) 



(6) 



is the field-free ground state of total energy Eq = 
Eq(Rq) + E VQ , where Eq(R ) denotes the electronic 
eigenenergy at nuclear equilibrium distance Rq, and E Vo 
is the vibrational eigenenergy of the nuclear Hamiltonian. 
The electronic part of this state is expressed as a Slater 
determinant of the orbitals (f>j(ri; Rq), while Xv {R) de- 
notes the vibrational ground state wave function, labeled 
by vibrational quantum number v§. Furthermore 
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— det [^(n; Ro)<h( r X R o) ■ ■ ■ $ q (r„, *)] 



,+ 



x X +(i?)expH£+t). 



(7) 



is the state of the molecule interacting with the laser field 
and E + = E e ' + (Ro) + E+ is the total energy of the resid- 
ual molecular ion. We also introduced <& q , which denotes 
the state of the active electron, and we neglect the inter- 
action of this electron with the residual molecular ion. 
Then the state of the active electron is given by a Volkov 
wave 



1 



(q+A(t))-r 



q • a sm(u>t) — — - sin(2wt) 



U p )t 



(8) 



with quiver radius a 
tential U p = -Ag/4. 



Aq/lu and ponderomotive po- 
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Using the Slater-Condon rules the expression in tron coordinate of the active electron 
Eq. JSJ reduces to an expression involving only the elec- 



do(t) = ~iJ2 I dt ' I d q(^ok) (Mr;RQ)\de\^(r,t))exp[i(E P {R Q )+E, - E e >+ (R ) ~ E+J t] 

x ( Vi \v ) ($ q (r',<')l V F (t') \Mr';Ro)) exp [-i {E%(R ) + E Vo - E e >+(R ) - E+) t'] . (9) 

I 



To simplify notation, we have dropped the index of the 
electron and consequently introduced d e — er, Vf(0 = 
F(i') • r'. Also we denote the HOMO wave function by 
cj) (r), and 




dRxt*{R)Xu {R) (10) 



denotes the Franck-Condon (FC) factor. 



A. The harmonics and the three-step picture 

The harmonic spectrum is calculated from the Fourier 
transform of the dipole moment do (t) , which implies that 
the intensity of the iVth harmonic is proportional to the 
norm square of the Nth Fourier component, djy- This 
Fourier component is investigated by extending the tech- 
nique of Ref. to the molecular case 



d N = ijf dteiq>{iNwt)d (t) = -J2J2f J o dt fJ Q dt ' 

x (27r) 2 K|^)(0 o (r)|exp( l ^)4|$ K , l (r,t))exp[z(i? o e (i? o ) + ^ o - £ e '+(i? ) - £+) t] 
x ( fih > ( (r',t')\V F (t')Mr') ) exp [i {E e >+(R Q ) + E+ - E e Q (R ) - E Vo ) t'] . (11) 

I 



An interpretation of Eq. (|llf) is possible in terms of 
the three-step picture of HHG. The matrix element 
( ^k"* (r', t') I Vfr (i') |<^o ( r ') ) describes the transition from 

k 

the HOMO and into the field-dressed Volkov state with 
momentum of magnitude 

K% = yj2(kw + ES{Ro) + E va - E*'+(Ro) - E+ - U p ). 

The direction of the momentum is parallel or anti-parallel 
(prescribed by the positive or negative sign of a de- 
fined below) to the polarization vector of the laser field. 
This momentum emerges after performing the q inte- 
gration in Eq. ©, provided that the laser field is suf- 
ficiently strong in the sense that the quiver radius, cxq, 
is large compared to the extent of the molecule. Fol- 
lowing |l5|. we include only terms corresponding to real 
values of because of the clear physical interpreta- 
tion of Eq. (|11|) . The exponential of the last line of 
Eq. I|llf> describes the energy balance in this initial step, 
and {vi\vo) is the FC factor associated with the transition 
in the nuclear degree of freedom. The factor l/Lo(t,t') 
with L(j{t,t') = <raa(smujt' — sin cut) and a — ±1 to as- 
sure Re(Lo) > comes from the Volkov propagator and 



describes the decrease of the amplitude of the contin- 
uum wave as the latter expands into three-dimensional 
space. Finally, the matrix element ( (/)o(r)|rf e |$ K "i (r, t) ) 

k 

with associated phase factors and FC factor describes the 
laser-assisted recombination (LAR) step of the three-step 
model. We note that, in our model, the expression of 
Eq. i|ll|) predicts a correlation of the momentum of the 
continuum electron with the nuclear motion of the resid- 
ual molecular ion, through the oscillating time-dependent 
exponential factors. 

B. The harmonic spectrum from arbitrarily 
oriented molecules 

In order to be able to evaluate Eq. I|ll|) for molecules 
of arbitrary orientation with respect to the laser field we 
use a HOMO wave function, where the angular part is 
expanded onto spherical harmonics. In the body-fixed 
frame, where the z axis is chosen along the internuclear 
axis, the HOMO wave function is then given by 

^(r)=^F,, m (r)^(f). (13) 

l 
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We also need the asymptotic form of this expression 

4/(t) ~ J2 C l,mr Z/K - 1 exp(-«r)y ; m (f). (14) 
i 



Here k = y/—2Eb, E^ is the binding energy of the active 
electron in the HF ground state wave function and Z is 
the net charge of the molecule, when the HOMO electron 
is removed. The t' integral of Eq. corresponds to the 
ATI step. It is well known that in this process the HOMO 
electron escapes to the laser-dressed continuum, when the 
distance to the nucleus is large ^^|. As a result we can 
employ the asymptotic form of the HOMO in the t' inte- 
grand in Eq. (|lip. and the integral will subsequently be 
evaluated using the saddle-point method. This evalua- 
tion proceeds as in Ref . |19| , and the HOMO as given by 
Eq. (|14|l is replaced in favor of the momentum space wave 
function. It is sufficient to know the latter near its sin- 
gular points, where it is given by the particularly simple 
expression Eq. (|A4(I . as described in Appendix A. The t 
integration occurring in Eq. corresponds to the LAR 
step and is performed numerically. Previous publications 
show that in the LAR step it is inadequate to approxi- 



mate the HOMO with the asymptotic form, because the 
recombination takes place at small separations of the ex- 
cited electron and the molecular core HH3. Also here 
we prefer to replace the space wave function by from 
Eq. (|13[) by the momentum space wave function, which 
in the body hxed frame may be written in the following 
way 

^(q) = ^G ; , m ( g )r ; m (q), (15) 

where the G; im 's are obtained numerically. Because the 
HOMO wave functions are expressed in terms of spher- 
ical harmonics, we can easily obtain the wave functions 
of the molecule in the laboratory frame, where the in- 
ternuclear axis does not necessarily coincide with the z 
axis. For a molecule, where the internuclear axis is ro- 
tated with respect to the z axis of the laboratory frame, 
the rotation is characterized by the three Euler angles 
= (a, (3, 7) [21]], and we arrive at the following expres- 
sion of the iVth harmonic for a molecule of orientation S% 
with respect to the laser polarization 



in 



E E d. 



?2* 

m' 2 ,7712 



Vi k sp(k) 



(16) 



with 



T 



K 1 m 



q'=K fc i +A(t; 



(17) 



and 



(2tt) 2 f T exp[i(Nwt-S(t))] 



L o(t>tsp(.k)) 



q=K^*+A(t) 



(18) 



Here D l m , m . i = 1, 2 denotes the matrix elements of 
the rotation operator [2lJ and 

S(t) = kojt + K% ■ a Q sin(wt) + sin(2wt) (19) 

is the quasi-classical action. The index sp(k) denotes the 
saddle-points. For each k the saddle-points t'(k) are 
defined by the condition S'(t' sp (k)) = 0, and we use the 
ones with < Re(t' sp (k)) < T along with lm(t' sp (k)) > 0. 
Note that Eq. (|17f) is evaluated in the complex t' plane 
and the factor (±l) Zl in Eq. 1)17(1 corresponds to the lim- 
its ±ik of the size q' of the electron momentum at the 
saddle-points. Letting the laser polarization define the z 
axis of the laboratory frame, we have both the momen- 
tum of the Volkov wave as well as the vector potential 



parallel to the z axis, so that Y l ^q') oc <5 m < i!0 in Eq l|17jl. 
Equation l(TBl) is evaluated along the real t axis, and when 
the polarization vector of the field is directed along the 
z axis e • V q reduces to a partial derivative with respect 
to the size q of the momentum. As mentioned earlier 
we consider only values of k, where is real, which 
serves to define a lower limit k m i n of the sum. The up- 
per limit k m ax is determined by requiring convergence 
of the calculated value of djv- Equation (fTfi|) along with 
the accompanying definitions of Eqs. (|17fl - (|19fl summa- 
rize the formulation of our model. A cruder estimate 
of the harmonic generation can be obtained by using the 
asymptotic form of the HOMO wave function in both the 
ATI and the LAR step. In this approximation the Nth 
harmonic is given by 
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with 



at T - ., ttt\ (— «e ■ Vq) 



r 



1 Tjh + 2 + Z/K) > 

3T -f; 2 w 



i«y 2'2+ 1 k 2 + z / k T(l 2 + |) 



_ ,'Z 2 + 2 + Z/ K ; 2 + 2 + Z/k 1 3 
x 2 Fi| , +2 ;/2+ 2 ; 



K / 



(21) 



q=K"*+A(t) 



where all functions are analytically known and where the 
Cim parameters are available in Table I (see also Table 
I in Ref. The asymptotic expression in Eq. 1|2U|) 

has the advantage that all geometric factors are pulled 
outside the integration. 

In Sec. IIVI we compare HHG spectra obtained with 
the two forms of Eqs. (|l(j|> and l|2(J|) . In that section we 
also show curves, which do not include nuclear motion. 
When the nuclei are clamped, the harmonic strengths are 
simply given by including only the first term in the sum 
over Vi in Eqs. I(16|) and 1)20(1 and replacing the Franck- 
Condon factors by unity. 



III. CALCULATIONAL DETAILS 

In this section we outline some calculational details to 
be used below. First we sketch the procedure used to ob- 
tain the Gi, m (q)'s and Ci, m 's of Eqs. JTHJ and (J2JJJ. The 
initial step in the procedure is to solve the HF equations 
fully numerically for the diatomic molecules [23 ■ In this 
manner the ground state orbitals are obtained, and are 
projected onto spherical harmonics to attain the radial 
functions F^ rn (r). The asymptotic angular coefficients 
are found by comparing the HF radial wave functions 
with the form Ci^ m r z l K ~ x exp(— nr), where the C; !m 's are 
used as fitting parameters. The G; jm (g)'s, on the other 
hand, are obtained from the Fourier-Bessel transform of 
the F,, m (r)'s. 

In the cases of N 2 and 2 we compare results of our 
model with experiments where the degree of molecular 
alignment was controlled by non-adiabatically creating a 
rotational wavepacket by a pump pulse and then moni- 
toring the harmonic yield as a function of the time delay 
to the intense probe pulse |lll Il2| . In order to perform 
a comparison with the experimental results, we need the 
time-dependent angular distribution p(t, /?) of the rota- 
tional wave packet created by the pump pulse. We refer 
to Ref. [2IJ for details about calculations of rotational 
wave packets. We assume that rotation is negligible un- 
der the probe pulse, and simply note at this point that 
once the angular distribution is available, the amplitude 



of the iVth harmonic at time t is given by 

d N (t)= [ d/3sm{P)p(t,(3)d N (P). (22) 
Jo 

where djv(/3) can be found using Eq. (|16|l for a probe 
laser parallel to the pump laser. We can then calculate 
the intensity of each harmonic from the rotational wave 
packet at every point of time. 

Table H] shows relevant information for the calculations 
presented in the next section. 



IV. RESULTS AND DISCUSSION 

A. The N2 molecule 

We begin by presenting a calculation of HHG from the 
N 2 molecule with a HOMO of a g symmetry. The laser is 
assumed to have an energy of u> = 0.057 a.u. correspond- 
ing to a wavelength of A = 800 nm, and the intensity is 
taken to be I = 2 x 10 14 W/ 'cm which is well below the 
saturation laser intensity [12|]. Such a choice of parame- 
ters is within the range of a TkSapphire laser. Figure ^ 
shows the relative strength of the harmonics as a func- 
tion of the harmonic order in a geometry, where the in- 
ternuclear axis of the molecule is aligned parallel to the 
polarization of the laser ((3 = 0°). Note that we refer to 
the norm square of the Fourier components (i.e. |g?at | 2 ) as 
the strength of the harmonics. There are three different 
curves on the figure, and the dotted curve corresponds to 
Eq. I|2U|) . i.e., the case where the asymptotic wave func- 
tion is used in both the ATI and LAR step of the HHG 
process, while the dashed and solid curves corresponds to 
Eq. (|16|l meaning that the improved wave function was 
employed in the LAR step. Calculations represented by 
the solid curve take nuclear motion into account. The 
harmonic spectrum of the N 2 molecule is different, when 
the asymptotic wave function is used everywhere as com- 
pared to the case, when the improved wave function is 
employed in the LAR step. This result shows that the 
quality of the HOMO wave function is of importance. It 
is also clear from Fig. ^ that the effect of nuclear mo- 
tion in the case of N 2 is negligible. This is because the 
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TABLE I: The molecular properties used in this work for evaluation of HHG. I p is the experimental adiabatic ionization 
potential, Ro is the equilibrium distance and B is the rotational constant of the molecule |24j|. The parallel and perpendicular 
polarizabilities an and a± of N2 and O2 are based on values from [2^1 and |2^|. respectively. We furthermore give the angular 
coefficients C; im from HF calculations. FC factors and vibrational energies can be found in the references indicated after each 
molecular species. 







Ip (eV) 


% (A) 


B (GHz) 


a„ (A 3 ) 


ax (A 3 ) 


Coo 


C2m 




N2 (<7 9 ) 


w 


15.58 


1.098 


59.647 


2.38 


1.45 


3.46 


1.64 


0.12 


o 2 K) 


[27] 


12.03 


1.208 


42.861 


2.3 


1.1 




1.04 


0.07 


H2 K) 


[28]° 


15.43 


0.741 








2.44 


0.14 




D2 (<T B ) 


f28] 


15.47 


0.742 








2.44 


0.14 





a Vibrational energies from Ref. |2'1 are used. 




Harmonic order N Orientation p (degrees) 



FIG. 1: (color online). Comparison of HHG from the N2 
molecule when using the asymptotic wave function in both the 
LAR and ATI step, and when the wave function used in the 
LAR step is improved. Only the solid curve includes nuclear 
motion. The signals have been normalized to the strength 
of the fifth harmonic obtained by using the improved wave 
function in the LAR step and including nuclear motion. The 
internuclear axis and the polarization are parallel (f3 = 0°). 



BO potential energy curves of the N2 molecule and the 
molecular ion are nearly parallel, which means that 
ionization of an N2 molecule in the ground state always 
leaves the residual molecular ion in the vibrational 
ground state. 

In Fig. [2] we study the orientation dependence of the 
HHG yield from N2. It is seen from the figure that har- 
monic generation for the parallel alignment (J3 = 0°) is 
more intense than the harmonic generation for the per- 
pendicular alignment ((3 = 90°), which reflects the orien- 
tation dependence of the ionization of the molecule (see 
Ref. |2(J and references therein). The detailed angular 
behavior depends on the harmonic order in considera- 



FIG. 2: (color online). Orientation dependence of the har- 
monics for the molecule Na- We have employed the improved 
wave function in the LAR step and included nuclear motion. 
All curves are normalized to their respective values at f3 = 0° . 



tion: a minimum is observed at approximately 40° for 
some orders. A similar behavior was recently predicted 
in Ref. 0, where the minima were located at a relative 
orientation close to 60°. 

Several experiments have measured the harmonic emis- 
sion from non-adiabatically aligned N2 molecules [ill IT^. 
E3. We can simulate these results as outlined in Sec. lllll 
To compare with the results of Ref. 

El 

we use the time- 
dependent angular distribution p{t, (3) corresponding to a 
rotational wave packet created by the experimental pump 
pulse, 60 fs, 800 nm, intensity 4 x 10 13 W/cm 2 , and the 
experimental initial rotational temperature of N2 of 30 
K. To model the results from Ref. JH| we obtain the 
time-dependent angular distribution with the experimen- 
tal pump pulse of duration 50 fs and with wavelength 800 
nm and intensity 6 x 10 13 W/cm 2 . In Eq. the value 
of djv(/3) is calculated using Eq. p6fl for a laser of wave- 
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Kanai et al. [11] 

?1 1 2 3 4 5 6 7 8 9 10 
Delay time (ps) 

FIG. 3: (color online). Comparison of theory and experiment 
in the case of N2. Panel (a) shows the average of the 21st 
to the 25th harmonics from a rotational wave packet of N2. 
Theoretical results are based on our model, while the results 
of Itatani et al. are experimental data from Ref. [T^ . In 
panel (b) the 23rd harmonic according to our model is com- 
pared with the experimental results of Ref. fTl| . Our results 
are normalized to the randomly aligned case. 

length A = 800 nm and intensity I = 2x 10 14 W/cm 2 , as 

111 Refs. mm 

We can then find the intensity of each 
harmonic from the rotational wave packet at every point 
of time. On Fig. [3] we compare model calculations with 
experimental results. The figure shows the evolution of 
the harmonic emission from the N2 molecule as a function 
of the delay time between the pump and probe pulse. In 
general the agreement between experiments and theory 
is very satisfactory. 

B. The O2 molecule 

We next consider the O2 molecule, which has a HOMO 
with TT g symmetry. The laser parameters are the same 




3 7 11 15 19 23 27 31 35 39 43 
Harmonic order N 



FIG. 4: (color online). HHG from the 2 molecule. The 
strength of the harmonics is normalized to the strength of 
the third harmonic obtained in the calculation, where an im- 
proved wave function is used in the LAR step and the nuclear 
motion is included. The internuclear axis is oriented at an 
angle (3 = 36° with respect to the internuclear axis. See text 
for further details. 

as in the preceding subsection. Figure 0] shows the rel- 
ative strength of the harmonics for /3 = 36°. As in the 
case of N2 it is important to use a better wave function 
than the asymptotic one in the LAR step. The effect 
of nuclear motion is more pronounced than for the N2 
molecule, because the Franck-Condon distribution of O2 
is broader than that of the N2 molecule (see Ref. 32] for 
more discussion on the implications on final state vibra- 
tional levels and possible control of these). 

We now turn to the angular dependence of the har- 
monics. The strength of some selected harmonics as a 
function of the orientation with respect to the laser field 
is shown on Fig. [SJ It is possible to understand the re- 
sults of Fig. from the symmetry of the HOMO. The ir g 
symmetry means that the HOMO has nodes along the 
internuclear axis as well as perpendicular to it, and this 
leads to a vanishing harmonic generation when (3 = 0° 
and (3 = 90°. Consequently the harmonics take on max- 
imal values for intermediate angels. All harmonics are 
most intense at orientations near 45°. It is also worth 
noticing that the precise location of the individual peak 
depends slightly on the order of the harmonic. This re- 
flects that higher harmonics probe regions closer to the 
core. 

For the O2 molecule, as well, several measurements of 
HHG from a non-adiabatically aligned sample have been 
carried out m> El> EH • We compare our model calcula- 
tions with the results presented in Ref. mi m the case 
of HHG from 02- In retrieving the time-dependent an- 
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FIG. 5: (color online). HHG from the O2 molecule as a func- 
tion of the orientation. The strength of the harmonics are 
calculated using Eq. I|16|l . All curves are normalized to their 
value at /3 — 45°. 



FIG. 6: (color online). Comparison of our model calculations 
with the experimental results from Ref. in the case of 
O2. We have calculated the intensity of the 23rd harmonic 
from a rotational wave packet of O2 normalized against the 
randomly aligned case. 



gular distribution we use the experimental pump pulse 
of 50 fs with wavelength 800 nm and intensity 6 x 10 13 
W/cm 2 . We calculate the value of d^iP) for use in 
Eq. (|22|l under assumptions corresponding to the exper- 
imental probe pulse with wavelength A = 800 nm and 
intensity I = 2 x 10 14 W/cm 2 . As in the case of N 2 
above, we normalize the harmonic intensity to the in- 
tensity of a randomly aligned sample in our theoretical 
results. The comparison is shown on Fig. and we ob- 
serve a fine agreement between theory and experiment. 



plies that a substantial part of the harmonic signal orig- 
inates from transitions, where the ionization potential is 
higher than in the fixed nuclei picture. As a result the 
whole process gets more unlikely and so the radiative 
emission is less intense. Very similar observations and 
conclusions hold true in the case of D2 . 



C. The isotopic molecules H2 and D2 

We now turn to a short discussion of calculations car- 
ried out for HHG from the isotopic molecules H2 and D2 . 
Harmonic spectra from these molecules provide a way to 
study the effect of nuclear motion, because the FC dis- 
tributions of both H2 and D2 are broad. Furthermore a 
comparison of the HHG signals from the isotopes con- 
tains information on how the nuclear masses affect the 
harmonic spectrum. 

In Fig. showing HHG from the H2 molecule, we have 
used the same laser parameters as in the case of N2 and 
2 studied above, i.e., 800 nm, 2 x 10 14 W/cm 2 . As seen 
from the figure, it is essential to include nuclear motion 
in the case of H2, and the general effect of the nuclear 
motion is to reduce the intensity of harmonic generation. 
This can easily be understood from the fact that the FC 
distributions of the isotope is broad. The broadness im- 



To see how the nuclear mass influences HHG, we com- 
pare the harmonic generation of the two isotopes for dif- 
ferent intensities at a laser energy of u> = 0.057 a.u. in 
a geometry, where the internuclear axes of the molecules 
are aligned perpendicular to the laser polarization. In 
order to be able to compare with results of Ref. @ 
we have furthermore included a calculation at inten- 
sity I = 4 x 10 14 W/cm 2 at energy uj = 0.0584 a.u., 
although our assumption of negligible depletion of the 
ground states is questionable in this situation. The re- 
sults are shown in Fig. |HJ In Ref. 8] a two-dimensional 
model was used and it was found that the heavier iso- 
tope leads to stronger harmonic generation. We can- 
not conclude something as simple from Fig. |H| Our non- 
conclusive picture of the relation between nuclear mass 
and HHG is supported by others , and is related to the 
strong correlation of the vibrational states of the residual 
molecular ion to the momentum of the excited electron 
as follows from Eq. (|llfl . 
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FIG. 7: (color online). HHG from the H2 molecule. Calcu- 
lations are performed under the setup mentioned in the text, 
and all curves are normalized to the strength of the fifth har- 
monic obtained by using the improved wave function in the 
LAR step and taking into account the effect of nuclear motion. 



V. SUMMARY AND CONCLUSION 

We have developed the strong-field approximation for 
HHG in diatomic molecules for calculating the harmonic 
signal from a single diatomic molecule of arbitrary orien- 
tation with respect to a monochromatic linearly polarized 
laser field. The model incorporates the effect of nuclear 
motion. 

We calculated the high-harmonic emission from the 
molecules N2, O2, H2 and D2. The two latter species 
were considered in order to study the effect of nuclear 
motion in detail. The effect of nuclear motion on these 
light molecules was significant and generally reduced the 
intensity of the harmonics. However, there seems to be 
no simple common relation between nuclear mass and the 
intensity of the harmonic generation. 

Turning to the N2, O2 cases, we successfully em- 
ployed the present model to reproduce experimental re- 
sults on HHG from non-adiabatically aligned N2 and O2 
molecules 

GUIII. The spectra for N 2 were very sensitive 
to the quality of the wave functions used to describe the 
molecule, while the effect of nuclear motion was small. 
The intensity of the harmonics was generally weaker, 
when the angle between the laser polarization and the 
internuclear axis was 90° as compared to the intensity at 
an angle of 0° . This orientation dependence characterizes 
diatomic molecules with a HOMO of a g symmetry. We 
also studied the orientation dependence of the harmonics 
for the O2 molecule, where the HOMO has ir g symmetry. 
As in the case of N2, the quality of the wave functions, 



FIG. 8: (color online). Comparison of harmonic generation 
for the isotopes H2 and D2 for different intensities at orien- 
tation f3 = 90°. The ordinate shows the strength of the ATth 
harmonic of D2 divided by the strength of the Nth. harmonic 
of H2. The laser parameters for the dash-dotted curve are the 
same as those used in Ref. 



used to model the O2 molecule was of importance, but 
furthermore we observed some effect of the nuclear mo- 
tion on the harmonic emission. The 7r g symmetry of the 
O2 molecule revealed itself through the angular depen- 
dence of the harmonics. In the case of O2 the harmonic 
generation was suppressed at orientations corresponding 
to both perpendicular as well as parallel alignment of the 
internuclear axis to the laser polarization. In support to 
earlier work 0, our results show that it is mainly the 
symmetry of the HOMO that determines the alignment 
dependence of the HHG. 

To conclude, we have presented a length gauge strong- 
field approximation for HHG in diatomic molecules. The 
theory is very versatile and readily adapted to a wide 
range of diatomic molecules, and straightforwardly gen- 
eralized to polyatomic molecules. The agreement with 
recent pump-probe experiments [rHll2| is excellent, and 
hence we have confidence in the predictive power of the 
model. Despite the good agreement, it is clear that the 
model poses important questions, in particular concern- 
ing the fact that our model is formulated in the length 
gauge. This seems to be advantageous for the initial ion- 
ization step in order to obtain the correct orientation 
dependence in the process, and it is also physically well- 
justified since the ionization is tunneling like and the 
escape into the continuum happens at large distances, 
where the length gauge interaction has large weight |16| . 
The physics in the recombination step, however, is in 
a sense completely the opposite. In this case an effec- 
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tive transfer of momentum has to take place between the 
recombining electron and the core. Such a process is ex- 
pected to be better described in the velocity gauge which 
naturally probes regions in space where the electron ex- 
periences abrupt changes (close to the nuclei). 

Another point is the question of how to treat the nuclei. 
The present formalism gives an attempt in that direction: 
within the Born-Oppenheimer approximation, we applied 
the Franck-Condon principle and observed significant ef- 
fects of the nuclear motion only for the lighter species H 2 , 
and D 2 . If it holds true then nuclear motion is essentially 
unimportant for all heavier systems, and it means that 
we may think of the nuclei as stationary fixed-in-space 
objects setting up an effective diffraction-like grating for 
the backscattered electron. Of course the out-come of 
the rescattering be it HHG or, e.g., future ionization, 
will be highly sensitive to this spatial arrangement, and 
this will undoubtedly lead to a higher degree of control 
of the harmonics. 
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APPENDIX A: THE ASYMPTOTIC 
MOLECULAR WAVE FUNCTION 

In this appendix we work out the asymptotic form of 
the momentum space wave function of the HOMO and 
describe, how this expression reduces to a well known 
function around singular points 0] . We start out in the 
body fixed frame of the molecule with the z axis chosen 
along the internuclear axis. The HOMO wave function in 
position representation must then follow the asymptotic 
Coulomb form 

«^(r) ~ ^G, m r z / K - 1 exp(- K r)r™(f), (Al) 
i 

where f is a unit vector in direction of r, — n 2 /2 = Ej, is 
the binding energy of the HOMO-electron, and Z is the 
net charge of the molecule, when the HOMO electron 
has been removed. ^ m (f) is a spherical harmonic. The 
summation is over the set of l's such that C;. m ^ 0. In 
this article we are concerned with terms corresponding 
to 1 = 0,2,4 (c.f. Tabled). 

The Fourier transform of Eq. IjAlf) yields the momen- 
tum space wave function 



„ [1 + 2 + Z/k 1 + 2 + Z/k 1 , 3 fq\ a \ 

x ^[ — r^' — 2-^ + 2 ; * + 2 ; -(J ) 

I 



where q is a unit vector in the direction of q. By 
evaluating Gauss' hyper geometric functions of type 
2 F 1 (a, a+ 1/2; I + 3/2; z 2 ) and then insert a = (I + 2 + 

Z/k)/2 and z 2 — ~ (q/n) 2 we can obtain an explicit rep- 
resentation of the Fourier transform in Eq. I|A2(I . To de- 



termine an expression of Gauss' hyper geometric function 
we have used formulas (15.1.9) and (15.1.10) in Ref. [34| 
along with a recursion relation coming from (9.137,1) in 
Ref. • Consequently, we obtained expressions of type 



2Fl (.,. + 1. * + _ ( _ ir (2n-l)(2„-3)...3 (1 - 2 T 2 ° 

x(2a- (2n- l))" 1 (2a- (2n - 2))~ 1 (2a - (2n - 3))" 1 . . . (2a - l)" 1 



[*(! + *)]' 



' 2 ("-^2 (1 + z)2 -« 



[z(l-z)\ 



(A3) 



where ^ (± , 2 (n — 1 ) ) is a polynomial of degree 2 (n — 1) . 



The first five of these polynomials are given in Table ITU 
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TABLE II: Polynomials of Eq. JX3j for use in the asymptotic form of the HOMO momentum space wave functions. 



^(±,2(n- 1)) 



l±(2a- l)z + (2a - 2)z 2 

3 ± 3(2a - l)z + (4a 2 - 7)z 2 ± (8a 2 - 18a + 7)z 3 + (4a 2 - 12a + 8)z 4 

15 ± 15(2a - l)z + 6(4a 2 - a - 8)z 2 ± 2(4a 3 + 12a 2 - 55a + 24)z 3 + (24a 3 - 72a 2 + 6a + 57)z 4 
±(24a 3 - 120a 2 + 168a - 57)z 5 + 8(a 3 - 6a 2 + 11a - 6)z 6 

105 ± (210a - 105)z + 15(12a 2 - 4a - 29)z 2 ± 5(16a 3 + 24a 2 - 190a + 87)z 3 + (16a 4 + 160a 3 - 760a 2 + 200a + 699)z 
±(64a 4 - 160a 3 - 640a 2 + 1750a - 699)z 5 + (96a 4 - 640a 3 + 1140a 2 - 140a - 561)z 6 
- 560a 3 



±(64a 4 



640a 2 + 1750a - 699)z 5 + (96a 4 - 640a 
+ 1640a 2 - 1810a + 561)z 7 + 16(a 4 - 10a 3 



1140a 2 

35a 2 - 50a + 24)z 



We are also interested in the limit q — > ±in corre- 
sponding to singularities of the asymptotic wave func- 
tion, because this expression is used for the HOMO wave 
function in Sec.[n] when the t' integration is evaluated by 



J 



the saddle-point method. The continuation of Eq. (|A2|) 
to the complex q plane is straight-forward [3fij . Around 
the points of singularity Eq. (|A2|) can be written as 



$'(q) 



1/2 



(2k) z / k T(1 + Z/k) 

( g 2 + k 2)1+Z/k ■ 



(A4) 



r 



where ±1 corresponds to q — > For the /-values used 
in this work Eq. (|A4|) can be verified directly by calcu- 
lating limits of the explicit representation of the relevant 



Gauss' hyper geometric functions from Eq. I|A3(1 and sub- 
stituting these expressions into the terms in Eq. ljA2fl . 
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